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Abstract 

In this paper, we consider queueing systems where the dynamics are non- stationary and 
state-dependent. For performance analysis of these systems, fluid and diffusion models have 
been typically used. Although they are proven to be asymptotically exact, their effectiveness as 
approximations in the non-asymptotic regime needs to be investigated. We find that existing 
fluid and diffusion approximations might be either inaccurate under simplifying assumptions 
or computationally intractable. To address this concern, this paper focuses on developing a 
methodology based on adjusting the fluid model so that it provides exact mean queue lengths. 
Further, we provide a computationally tractable algorithm that exploits Gaussian density in 
order to obtain performance measures of the system. We illustrate the accuracy of our algorithm 
using a wide variety of numerical experiments. 

1 Introduction 

There are several applications of systems where the dynamics are state-dependent including the 
repairman problem, retrial queues, chemical reactions, epidemic models, communication networks 
with state-dependent routing, call centers, etc. Even assuming Markovian properties, analysis of 
state-dependent systems is difficult. Therefore, typically fluid and diffusion approximations are 
used to obtain performance measures of these systems [H O El 13 El 110] • These fluid and diffusion 
models are obtained by utilizing Functional Strong Law of Large Numbers (FSLLN) and Functional 
Central Limit Theorem (FCLT) which are well summarized in ? ] and ? ]. Using these models, one 
can investigate the asymptotic behavior of the system state which can be a good approximation 
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to the original system under certain specific conditions (e.g. heavy traffic, large number of servers, 
etc). The previous studies mentioned above have a common feature in that they utilize FSLLN 
and FCLT. They, however, have different aims, scenarios, and assumptions. This paper specifically 
focuses on the fluid and diffusion models for the non-stationary (i.e. time-dependent or transient) 
and state-dependent exponential queueing systems similar to those in ? ] and (author?) [6]. 
In fact, some rate functions to describe the system dynamics are of the forms, min(-, •) and max(-, •) 
which are not differentiable everywhere. Notice that to apply the seminal weak convergence result 
in (author?) [5], we require differentiability of rate functions which is not satisfied in most non- 
stationary and state-dependent queueing systems. To extend the theory to non-smooth rate func- 
tions, (author?) [6] proves the weak convergence by introducing a new derivative called "scalable 
Lipschtz derivative" and provides models for several queueing systems such as Jackson networks, 
multiserver queues with abandonments and retrials, multiclass priority preemptive queues, etc. In 
addition, several sets of differential equations are also provided to obtain mean values and the 
covariance matrix of the limit process. It, however, turns out that the resulting sets of differential 
equations are not computationally tractable in general and hence the theorems cannot directly be 
applied to obtain performance measures numerically. In a follow-on paper, (author?) [7] provides 
numerical results for queue lengths and waiting times in multiserver queues with abandonments and 
retrials by including an additional assumption to deal with computational tractability. Specifically, 
that paper assumes measure zero at a set of time points where the fluid model hits non-differentiable 
points, which enables applying Kurtz's diffusion models. However, as pointed out in (author?) 
j7], if the system stays in a critically loaded state (non-differentiable points) for a long time, also 
called "lingering", their approach may cause significant inaccuracy. Before describing our goal, it 
is worthwhile to summarize the above results and point out possible problems: 

• On one hand, (author?) [6j provides rigorous theory to obtain the fluid and diffusion models 
for the system having non-smooth rate functions. On the other hand, it is not possible to 
solve the resulting set of differential equations to obtain performance measures numerically. 

• The additional assumption of measure zero in (author?) [7] (also see ? ]) provides a com- 
putationally tractable way to obtain performance measures. However, when this assumption 
is not valid, it might cause inaccuracy in obtaining the results. 
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The goal of this paper is to devise a technique that strikes a balance between accuracy and com- 
putational tractability leveraging upon the fluid and diffusion models in (author?) [5]. To do so, 
we explain when the fluid approximations might fail to achieve accuracy from a different point of 
view than those considered in (author?) [6, 7]. Our method is irrelavent to the smoothness of rate 
functions and provides a condition to obtain exact estimation of mean values of the system state. 
We apply our methodology to several queueing systems including not only multiserver queues with 
abandonments and retrials considered in (author?) [6l [7] but also more complex queueing sys- 
tems such as multiclass priority preemptive queues (considered but not numerically investigated in 
(author?) [6j) and peer-based networks in multimedia distribution. Here, we emphasize that this 
paper does NOT aim at proving another weak convergence to a limit process but pursues providing 
a practically effective methodology to increase accuracy in performance measures such as mean 
values and the covariance matrix of the system state. Our paper has discriminating features from 
previous research in that we: 

• address possible inaccuracy of the fluid model which might occur irrespective of the smooth- 
ness of rate functions, 

• solve the fluid model directly by providing a methodology to estimate mean values exactly 
unlike previous research where the fluid model is unchanged and is complemented by the 
expected value of the diffusion model, and 

• devise an effective algorithm transforming the fluid and diffusion models, which achieves not 
only increased accuracy but also computational feasibility. 

We now describe the organization of this paper. In Section [2] we explain the fluid and diffusion 
models in (author?) [5j and (author?) [6 , and describe their limited applicability in practice. In 
Section [3j we provide a methodology to estimate exact mean values of the system state. However, 
this would not immediately result in a computationally feasible approach. For that, in Section |4j 
we explain our algorithm based on Gaussian density to achieve computational tractability and the 
benefits of using Gaussian density. In Section [5] we show how our proposed method works for the 
queueing system described in (author?) [7^ by comparing our method with theirs. In Section |6j 
we provide some numerical results for more complex queueing systems where we have reasons to 
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believe our methodology may not be accurate, and for these cases, we determine the performance of 
our approach. Finally, in Section [Tj we make concluding remarks and explain directions for future 
work. 



2 Recapitulating fluid and diffusion approximations 

Before explaining our results, we recapitulate fluid and diffusion approximations developed by 
(author?) [5] that we would leverage on for our methodology. As a matter of fact, the diffusion 
model developed in (author?) [5 is not directly applicable in many queueing systems because 
it requires differentiability of rate functions which is sometimes not satisfied. Therefore, we also 
briefly mention the result in (author?) [6] which extends the Kurtz's result to the model involving 
non-smooth rate functions. Further, it is worthwhile to note that for n G N, the state of the 
queueing system Xn{t) includes jumps but the limit process is continuous. Therefore, the weak 
convergence result that is presented is with respect to uniform topology in D (? ] and ? ]). 
Let X{t) be a d-dimensional stochastic process which is the solution to the following integral 
equation: 



{1,2,. ..,k} are constants, and /j's are continuous functions such that \fi{t,x)\ < Cj(l + for 
some Ci < oo, t < T and T < oo. Note that we just consider a finite number of /j's to simplify 
proofs, which is reasonable for real world applications. It is usually not tractable to solve the 
integral equation ([T]). Therefore, to approximate the X{t) process, define a sequence of stochastic 
processes which satisfy the following integral equation: 



Typically the process Xn{t) (usually called a scaled process) is obtained by taking n times faster 
rates of events and 1/n of the increment of the system state. This type of setting is used in the 




(1) 



where xq = ^(0) is a constant, 1^'s are independent rate 1 Poisson processes, h G Z'^ for i G 




^-^ n 

i=l 
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literature and is denoted as "uniform acceleration" in ? ], and (author?) [61 [7]. Then, the following 
theorem provides the fluid model to which {Xn{t)} converges almost surely as n ^ oo. Define 



F{t,x) = Y,hMt,x). 



(2) 



i=l 



Theorem 1 (Fluid model, (author?) [5|). If there is a constant M < oo such that \F{t,x) — 
F'{tjy)\ ^ M\x — y\ for all t < T and T < cxd. Then, lim„^oo ^n(i) = ^(t) <i-S- where X(t) is the 
solution to the following integral equation: 



X{t)=X0 + y2li I f^{s,X{s))ds. 

*=1 -^0 



Note that X{t) is a deterministic time- varying quantity. We will subsequently connect X{t) and 
X{t) defined in equation ([T]), but before that we provide the following result. Once we have the fluid 
model, we can obtain the diffusion model from the scaled centered process {Dn{t)). Define Dn{t) 
to be -y/n {Xn{t) - X{t)). Then, the limit process of Dn{t) is provided by the following theorem. 



Theorem 2 (Diffusion model, (author?) [5]). If fi's and F, for some M < oo, satisfy 



\fi{t,x)- fi{t,y)\<M\x-y\ and 



_d_ 

dx, 



F{t,x) 



<M, fori£{l,...,k} andO<t<T, 



then lim„_+oo Dn{t) = D{t) where D{t) is the solution to 

D{t) = y^k [' Jf,{s,X{s))dWi{s)+ [' dF{s,X{s))D{s)ds, 
Jo Jo 

Wi{-)'s are independent standard Brownian motions, and dF{t,x) is the gradient matrix of F{t,x) 
with respect to x. 

Remark 1. Theorem^ requires that F{-,-) has a continuous gradient matrix. Therefore, if we 
don't have such an F, then we cannot apply Theorem directly to obtain the diffusion model. 

Remark 2. According to (author?) f^, if D(0) is a constant or a Gaussian random vector, then 
D{t) is a Gaussian process. 
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Now, we have the fluid and diffusion models for Xn{t). Therefore, for a large n, Xn{t) is approxi- 
mated by 

X^{t) « X{t) + 



n 



If we follow this approximation, we can also approximate the mean and covariance matrix of Xn{t) 
denoted by and C ov\^Xn{t) , Xn{t)\ respectively as 

E\Dit)\ 

E[Xn{t)] « X{t) + L ^ (3) 



n 



r . s , ST Cov\Dit),D{t)\ 

Cov[Xn{t),Xn{t)] « ^ 4 

In equations ([s]) and Q, only is known. Therefore, in order to get approximated values of 

E[Xn{t)\ and Cov[Xn{t), Xn{t)\, we need to obtain E[D{t)\ and Cov[D{t),D{t)\. The following 
theorem provides a methodology to obtain and Cov\D{t),D{t)\. 

Theorem 3 (Mean and covariance matrix of linear stochastic systems, (author?) [1]). Let Y{t) 
he the solution to the following linear stochastic differential equation. 

dY{t) = A{t)Y{t)dt + B{t)dW{t), y(0) = 0, 

where A{t) is a dxd matrix, B{t) is a dxk matrix, and W(t) is a k-dimensional standard Brownian 
motion. Let M{t) = E\Y{t)\ and T.{t) = Cov[Y{t),Y{t)]. Then, M{t) and S(f) are the solution 
to the following ordinary differential equations: 

^M{t) = A{t)M(t) 

jT.{t) = A{t)T.{t) + T.{t)A{ty + B{t)B{t)'. (5) 
Corollary 1. //M(0) = 0, then E[M{t)] =0 fort>0. 

By Corollary [l| if D{0) = 0, then E[D{t)] = for f > 0. Therefore, if X{0) = X{0) = xq, then we 



6 



can rewrite ([s]) to be 



E[Xn{t)\ « X{t). 

Recalling Remark[T] the diffusion model in (author?) [5i requires differentiability of rate functions. 
Otherwise, we cannot apply Theorem |2j To address this problem, (author?) [6] introduces a new 
derivative called "scalable Lipschitz derivative" and proves the weak convergence using it. Unlike 
the result in (author?) [5^, it turns out that the diffusion limit may not be a Gaussian process 
when rate functions are not differentiable everywhere. In (author?) [^, expected values of the 
diffusion model may not be zero (compare it with Corollary [T]) and could adjust the inaccuracy in 
the fluid model (see (author?) [7 ). The resulting differential equations for the diffusion model, 
however, are computationally intractable. For example, in (author?) [6], one of the differential 
equations has the following form: 

rendering it to be intractable. 

Therefore, (author?) [7 , as we understood, resorts to the method in (author?) |5] by assuming 
measure zero at non-smooth points to avoid computational difficulty. As described in Section [T| in 
this paper, our objective is to give the fluid and diffusion models a fresh look from an alternative 
perspective, and suitably adjust them for non-asymptotic scenarios. This is presented in the next 
section. 



3 Adjusted fluid model 

In this section, we first explain the possibility 
system state using the fluid model. Then, we 
mean values. For those, we consider the actual 
by the following theorem. 



of inaccuracy when obtaining mean values of the 
provide an adjusted fluid model to estimate exact 
integral equation to get the exact value of £^[X(t)] 
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Theorem 4 (Expected value of X{t)). Consider X{t) defined in equation |Ip. Then, for t < T, 
E[X{t)~\ is the solution to the following integral equation. 



E[X{t)]=xo + Y,k E\fi{s,X[s)) 



ds 



(7) 



Proof. Take expectation on both sides of equation (nl. Then, 



E[X{t)\ = xo + ^hE 



1=1 

k 



Yi(^l^ fi{s,X{s))ds 
X0 + / fi{s,X{s))ds 

i=i Lio 



since ^(O's are nonhomogeneous Poisson processes 







i=l 



xo + ^^^j / E /j(s,X(s)) by Fubini theorem in (author?) ^3]. 



Therefore, we prove the theorem. 



□ 



Comparing Theorems [l] and |4j notice that we cannot conclude that X{t) in Theoremjljand ii^[X(t)] 
in Theorem |4] are close enough since E[fi{t, X{t))~\ / fi(^t, E[X{t)]y In some applications, /j's 
might be constants or linear combinations of components of X{t). In those cases. Theorem |4] and 
the following corollary imply that the fluid model would be the exact estimation of mean values of 
the system state. 

Corollary 2. If fi{t,x)'s are constants or linear combinations of the components of x, Then, 

E[X{t)]=X{t), 



where X{t) is the solution to |jp and X{t) is the deterministic fluid model from theorem^ 

Proof. Using linearity of expectation in ? ], we can obtain the same integral equation for both 
E[X{t)] and X{t). □ 

However, if we have different forms of fiS where E[fi{t, X{t))\ ^ fi{t, E[X{t)]) , then the fluid 
model would be inaccurate. As seen in Section [2] the fluid model does not require differentiability 
of rate functions in both (author?) [5j and (author?) [6 . In (author?) [6], the diffusion model 
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can contribute to mean values of the system state. However, as seen in equation ([6|, the differential 
equations to obtain mean values of the diffusion limit are not computationally tractable. Even if 
they are numerically solvable, mean values of the diffusion limit is zero by the time the fluid limit 
hits a non-differentiable point for the first time. We will show in Section |6] that inaccuracy begins 
to occur before the fluid limit hits that point. Therefore, we approach this problem in a different 
point of view. 

The basic idea of our approach is to construct a new process {Z(t)) so that its fluid model is exactly 
same as mean values of the original process X{t) as described in Theorem |4j Define a set F of all 
distribution functions that have a finite mean and covariance matrix in R*^. This set is valid for 
the fluid model since conditions on /j's guarantee that < oo and |Cow[X(t), | < oo 

for all t < T. Define a subset Fq of F such that any / G Fq has zero mean. We call an element of 
Fq a "base distribution" for the remainder of this paper. 

Proposition 1. E[fi{t, X{t))] can be represented as a function of E[X{t)] for t <T. 

Proof For fixed to ^ T, suppose the distribution of X^to) is F. Then, F €¥. For F S F, we can 
always find Fq G Fq such that F{x) = Fq{x — /i) where = E[X(to)] = j-^d xdF. Then, 

E[fi{ta,X{ta))] = [ fi{to,x)dF 

= / fi{to,x + fi)dFo. 

Since the integration removes x, by making to and /x variables (i.e. substitute to and /i with t and 
/i(t) respectively), we have 

E[fi{t,X{t))] = gi{t, fi{t)), for some function gj. 

□ 

Remark 3. Proposition^ does not mean that fi{-) completely identifies the function gi{-,-). In 
fact, the function gi{-,-) might be unknown unless the base distribution is identified but we can say 
that such a function gi{-,-) exists. 
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For t < T, let = E[X{t)]. Let gi{t,fi{t)) = E[fi{t, X{t))] for i E {l,...,k}. Then, we can 
construct a new stochastic process Z{t) which is the solution to the following integral equation: 



Z{t) = zo + J2liYi(^j\i{s,Z{s))dsy 



(8) 

Based on equation (|8]), define a sequence of stochastic processes {Zn{t)} satisfying 

(9) 



Zn{t) = xq + Y^^hYi(^j ngi{s,Zn{s))ds 



Next, we would like to obtain the fluid model for Zn{t). Before doing that, we, however, need 
to check whether the functions giS satisfy the conditions to apply Theorem [T] Following lemmas 
provide the proofs that gis meet the conditions. 

Lemma 1. // \ fi(t, x)\ < Cj(l + |x|) for t <T, then gi{t, x) 's satisfy 

\9i{t,x)\ < Di{l + \x\) for some Di < CO. 

Proof. To prove this lemma, we need to show that < + |^ for K < co and 

t <T. We first show it in the one-dimensional case and then extend it to the d-dimensional case. 
Let, for fixed to <T, X = X{to) having mean fj, and variance o"^, and fi{X) = /j(to,-^(to))- Then, 
by Cauchy-Schwarz inequality, 

E[\X\]<^/^X^=^/ll?T^<\^i\+CJ<D{l + \^l\) for il) = max(l, cj). (10) 

Now, we have the one-dimensional case and can move to the d-dimensional case. Suppose X has a 
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mean vector fj, and a covariance matrix S such that X = (xi, . . . , Xd)', /U = (/^i, • • • , fJ-d)' ■ Then, 



E[\X\] = E 



< E 



i=l 



=Y,E[\x,\] 



< D(d + ^\fii\] by ((To]) for L> = max(l,f7i,...,crd) 
^ i=i ' 



< D\^d + d 
= Dd{l + \fi\). 



A I by Cauchy-Schwarz inequahty 

\ i=i 



(11) 



Now we have < Kyi + for the d-dimensional random vector X where K = Dd. 

Then, 



E[fiiX)] < E\\fi{X)\\ <a + CiE[\X\] from assumption 



< Ci + QK{1 + \fi\) < A(l + for Di = Ci + CiK by equation ([TTJ 



Note gi{tQ,^) = Since |E| is bounded on t < T, if we make to > arbitrary, we prove 

the lemma. □ 



For the next lemma, we would like to define 



G{t,x) = ^kgiit,: 



1=1 



Lemma 2. For t <T, if \ fi{t,x) — fi{t,y)\ < M\x — y\, then gi{t,x) 's satisfy 



\gi{t,x) - gi{t,y)\ < M\x - y\, 



and if\F{t,x) — F{t,y)\ < M\x — y\, then G{t,x) satisfies 



\G{t,x)-G{t,y)\<M\x-y\. 



(12) 



Proof. For fixed to < T, let X = X(to) and Y = Y^to) and suppose X and Y have a same base 
distribution Hq (we use H instead of F to avoid confusion with F in ([2])) where E[X] = /ii and 
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E\Y] = ^2- Then, the distribution Hi of X and H2 of Y satisfy 



Hi{x) = Hq[x — ^i), and 
H2{y) = Ho{y-^i2), 



respectively. Now, we have 



E[F{X)] -E[F(Y)] 



f F{x)dHi- [ F{y)dH2 



By transforming variables. 



E[F{X)] - E[F{Y)] = [ F{x + fii)dHo- [ F{y + f,2)dHo 



[ {F{x + fii)-F{x + fi2))dHo 



by linearity. 



< 



dHn 



< M \fj.i — fi2\dHo = M\fii — fi2\ by assumption. 

jRd 

Note G{to, /ui) = E[F{X)] and G{to, H2) = ^[^C^^)] • Then, by making to > arbitrary, we prove 
the second part, i.e. if \F{t,x) - F{t,y)\ < M\x - y\ then \G{t,x) - G{t,y)\ < M\x - y\. We can 
prove the first part, i.e. if \ fi{t, x) — fi{t, y)\ < M\x — y\, then \gi{t, x) — gi{t, y)\ < M\x — y|, in a 
similar fashion and hence we have the lemma. □ 

Lemmas [1] and [2] show that if /j's satisfy the conditions to obtain the fluid limit of then giS 

are also eligible for the fluid model of Zn{t). Therefore, we are now able to provide the adjusted 
fluid model based on Lemmas [T] and |2j 

Theorem 5 (Adjusted fluid model). Assume 



\h{t,x)\ < C,(l + |x|) foriG{l,...,k}, 
\F{t,x)-F{t,y)\ < M\x-y\. 



(13) 
(14) 
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Then, lim„^oo Zn{t) = Z{t) a.s., where Z{t) is the solution to the following integral equation: 



Z{t)=XQ + ^lil gi{s,Z{s)) 
^=l 



ds. 



and furthermore 



(15) 



Z{t)=E[X{t)]=xo + Y.^^ / E\fi{s,X{s)) 



i=l 



ds. 



Proof. From Lemmas [T] and [2| (13) and (14) imply 



(16) 



\gi{t,x)\<Di{l + \x\) and \G{t,x) - G{t,y)\ < M\x - y\ 



Therefore, by Theorem [T] we have equation (15), and by definition of gi{t,x)^s, we have equation 
(161). □ 



In Theorem [s] we have the same conditions for /j's and gi^s, and gi{t,x)^s do not guarantee 



E[gi{t,X{t))] = gi{t,E[X{t)]) either. However, comparing equation ( |l6[ ) with equation 
Theorem |4| we notice that Theorem [5] via equation ( 16 ) could provide the exact estimation of 
E[X{t)]. 

Though Theorem [5] provides the exact estimation of £'[X(t)] , we should identify the functions gi^s 
in order to obtain these values numerically or analytically. The (/j's, however, cannot be identified 
unless the base distribution is known, which forces us to develop an algorithm to find gi^s. The 
following section will describe our Gaussian-based method which would also be useful to adjust the 
diffusion model. 



4 Gaussian-based adjustment 

In Theorem [5} we encounter a fundamental problem in finding gi^s, i.e. we need to characterize 
the distribution of X{t). There, however, is no clear way to find the exact distribution of X{t) 
in general. Therefore, our proposed method starts with assuming the distribution of X(t). Recall 
that in Section [2| X„(t) is approximated by a Gaussian process when xq is a constant for a large 
n. Though it is not true for n = 1, we use a Gaussian density function to obtain gi^s with zq = xq 

13 



since using the Gaussian density function provides following three benefits: 

1. In many applications, empirical densities are close to Gaussian density even if rate functions 
are not differentiable (see (author?) [8j, (author?) [7J. 

2. Gaussian distribution can be completely characterized by the mean and covariance matrix 
which can be obtained from the fluid and diffusion models. 

3. By using Gaussian density, g^'s can achieve smoothness even if /j's are not smooth, which 
enable us to apply Theorem [2] directly. 

The third benefit is not obvious and hence we provide the proof of that. 

Lemma 3. Let gi 's be the rate functions of Z{t) obtained from Gaussian density. Then, gi 's are 
differentiable everywhere. 

Proof. Define 




Using Gaussian density. 




For j G {1, . . . ,d}, since (f>{x,y) is differentiable with respect to Xj and \fi{t,y)-^(l){x,y)\ is inte- 
grable. 



d 

dxj 





fi{t-:V)~i — 4'[x,y)dy by applying Theorem 2.27 in (author?) [7, (17) 



where Xj is j component of x. 

Therefore, gi is differentiable with respect to Xj. 



□ 
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Now, we have gi{-,-)^s which are differentiable. Then, we can apply Theorem |2] to obtain the 
diffusion model for Zn{t). 

Proposition 2 (Adjusted diffusion model). Let gi{-, ■) 's be the rate functions in Z{t) obtained from 
Gaussian density. Define a sequence of scaled centered processes {Vn{t)} for t <T to be 



Vn{t)=Vn{Zn{t)-Z{t)), 



where Zn(t) and Z(t) are solutions to equations Im and (15) respectively. If fi{t,x)'s and F{t,x 



satisfy equations (13) and (14) respectively, then lim„_^oo 

Vn{t) = V{t), where 



V{t) = y^h [' Jgi{s,Z{s))dW,{s)+ f dG{s,Z{s))ds, 

Jo 

Wi{-)'s are independent standard Brownian motions, and dG(t, Z{t)) is the gradient matrix of 
G{t,Z{t)^ with respect to Z{t). Furthermore, V{t) is a Gaussian process. 



Proof. From definition of G{t, x) in ( |12[ ), we can easily verify that G{t, x) is differentiable by Lemma 
[3] and hence \G{t,x) - G{t,y)\ < M\x - y\ implies 



d_ 
dx, 



G{t,x] 



< Mi for some Mi < oo, t <T, and i £ {1, . . . , d}. 



Therefore, by Theorem [2] we prove this proposition. 



□ 



Corollary 3. If fi's are constants or linear combinations of the components of X{t). Then, 



X{t) = Z{t) in distribution. 



Proof. Using the linearity of expectation, we can verify gi{t, x) = x) for iG{l,...,A;}. □ 

Now we have the adjusted fluid and diffusion models by utilizing Gaussian density. Therefore, 
instead of assuming measure zero at a set of non-differentiable points (as done in (author?) [7]), 
we compare the adjusted models with the empirical mean and covariance matrix. Note when 
we explain Theorem^ we do not consider S(t), the covariance matrix of X{t). However, from 
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Gaussian density, we know that characterizes the base distribution and it can be obtained 
from Proposition [2] Therefore, we rewrite y^'s to be functions of t, Z{t), and i.e. 

gi{t,Z{t)) ^ gi{t,Z{t),^it)) for i G {!,..., A;} and (18) 
G{t,Zit)) ^ (19) 

Proposition 3 (Mean and covariance matrix). Let Y{t) = Z{t) + V{t). Then, 

E{Y{t)) = Z(t) and (20) 
Cov{Y{t),Y{t)) = Cov{V{t),V{t)) =T.{t). (21) 

The quantities Z{t) and S(t) are obtained by solving the following simultaneous ordinary differential 
equations with initial values given by Z(0) = xq and S(0) = 0.' 

d - ^ 

-z{t) = Y.im{t,z{t),nt)), (22) 

i=l 

^m) = A{tMt) + ^it)Aity + Bit)Bity, (23) 



where A{t) is the gradient matrix of G(t, Z{t),T,(t)^ with respect to Z{t), and B{t) is the dx k 
matrix such that its i^^ column is li\J giit, Z{t), S(t)) . 

Proof. Since V^(0) = 0, from Corollary [T| we have (20) and ( pTj ). By rewriting (15) in Theorem [5] 



as a differential equation form, we have l\22h, and by Theorem 3 we have (23). Note that since 



both Z{t) and S(t) are variables, we should solve (22) and ([23| simultaneously. □ 

In conclusion, we define an adjusted process Z{t) in Section [s] to obtain the exact £'[X(t)] for X{t) 
process which is the state of a non-stationary and state-dependent queueing system. It, however, 
is not possible to obtain such gi^s and hence in this section, we provided an algorithm by utilizing 
Gaussian density. From this, the limit process turns out to be a Gaussian process. We recognize 
that this is not true for the original process. As mentioned in Section [T| however, this paper does 
not pursue finding the exact distribution of the original process but proposes an effective way to 
estimate mean values and the covariance matrix of the original process. Therefore, in the following 
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Figure 1: Multiserver queue with abandonment and retrials, (author?) [7] 

sections, by means of numerical examples, we illustrate our methodology and show its effectiveness. 



5 Multiserver queues with abandonments and retrials 

Multiserver queues with abandonments and retrials are extensively studied in the literature since 
they are used to model an important application, namely "call centers" (e.g. [? ? El? , EO]). In 
this section, therefore, we provide in-depth explanation of how our approach works in this queueing 
system by numerical examples. Figure [T] illustrates a multiserver queue with abandonments and 
retrials described in (author?) [6l [7j. There are nt number of servers in the service node at time 
t. Customers arrive to the service node according to a nonhomogeneous Poisson process at rate A^. 
The service time of each customer follows a distribution having a memoryless property at rate . 
Customers in the queue are served under the FCFS policy and the abandonment rate of customers 
is Pt with exponentially distributed time to abandon. Abandoning customers leave the system with 
probability pt or go to the retrial queue with probability 1 — pt. The retrial queue is an infinite- 
server-queue and hence each customer in the retrial queue waits there for a random amount of time 
with mean l/fJ-f and returns to the service node. 
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Let X{t) = {xi{t), X2{t)^ be the system state where xi{t) is the number of customers in the service 
node and 2:2 (t) is the number of customers in the retrial queue. Then, X{t) is the unique solution 
to the following integral equation: 



xi{t) = xi(0) + Xsds^ ^"^^(/o An,)/ija!s) 



Then, following the notation in Section [2| we have, for x = {xi,X2) and t <T, 

fi{t,x) = Xt, 

f2{t,x) = 1^X2, 

h{t,x) = filixiAnt), 

fi{t,x) = f3t{l-pt){xi-nt)^, and 

/5(t,x) = I3tpt{xi - nt)^ . 

We can verify that all /j's satisfy the conditions to apply Theorem [TJ However, we cannot apply 
Theorem [2] directly since /s, /4, and /s are not differentiable at xi = n^. To resolve this, (author?) 
[7] assumes measure zero at a set of time points when the fluid limit hits the non-differentiable 
points and apply Theorem [2| (author?) addresses that for a system of a fixed size, assuming 
measure zero works well when xi{t) does not stay too long near the critically loaded phase. It 
also provides the actual form of differential equations for the diffusion model in (author?) [6] 
which, in fact, are not computationally tractable, e.g. see equations (4.1) and (4.2) in [7 . As 
mentioned in Section |3] we approch the problem under a different point of view. Notice that in 
addition to their non- differentiability, /s, /4, and /s do not satisfy E[fi[t,X{t))] = fi{t,E[X{t)\) 
either. Therefore, we would like to apply Theorem |5] to obtain £'[X(t)] exactly. Recalling Section 
|4j however, obtaining exact gis is not possible and hence we obtain gis from Gaussian density as 
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follows: 



gi{t,x) = Xt, 

92{t,x) = f4x2, 

93{t,x) = f4{nt + {xi-nt)<^{nt,xi,ai,) -al^(/){nt,xi,ai,)), 

gA{t,x) = I3t{l-pt){{xi-nt){l-^{nt,xi,ai^)) +al^(t){nt,xi,aiS^, and 

g^{t,x) = I3tpt{{xi-nt){l-^{nt,xi,ai^)) +al^4){nt,xi,(JiS^, 

where $(a, 6, c) and (/)(a, 6, c) are function values at point a of the Gaussian CDF and PDF respec- 
tively with mean h and standard deviation c. 

Since fi{t,x) and f2{t,x) are constant and linear with respect to x respectively, gi{t,x) = fi{t,x) 
and g2{t,x) = f2{t,x). The derivation of other gi{-, •)'s is straightforward but requires some com- 
putational efforts and hence we provide the details in Appendix [A} Note 53, 54, and 55 include 
(Tij which is currently treated as a function of t but is used by the adjusted diffusion model (see 



equations (18) and (19)). 

For a fixed size of the system, both our proposed method and the method assuming measure zero 
do not guarantee the exact estimation of the system state. However, these two methods provide 
computational tractability. Therefore, we compare our method against the method assuming zero 
in (author?) [7 . We conduct simulations under the similar settings in (author?) [7 . We use 
5,000 independent simulation runs and compare the simulation result with both our method and 
the method assuming measure zero. We use the constant rates for the parameters except the arrival 
rate. The arrival rate alternates between 45 and 55 every two time units. Figures |2] and |3] show the 
estimation of mean values from one experiment. The number of servers {nt) is 50 and the service 
rate of each server is 1. As seen in Figure [2] the number of customers in service node (xi(t)) stays 
near the critically loaded point for a long time. As (author?) [7 points out, the method assuming 
measure zero shows significant difference in estimating i?[x2(t)] . On the other hand, our proposed 
method provides accurate results. When we see the estimation of the covariance matrix, we notice 
the similar results as the estimation of mean values. As seen in Figure [3j the method assuming 
measure zero causes "spikes" as also pointed out in (author?) [7]. Our proposed method, however. 
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n =50,H|=1 .|ij^=0.2, p=0.5, =45, i.f=55 n =50,H| =1 .11^=0.2. p =0.5, =45, ?i|^=55 




Time Time 



(a) Mean numbers by assuming measure zero (b) Mean numbers by our proposed method 

Figure 2: Comparison of mean values, E[X{t)] 
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provides reasonable accuracy and no spikes at all. To verify the effectiveness of our method, we 
conduct several experiments with different parameter combinations. Table [T] describes the setting 

Table 1: Experiments setting 



exp # 


svrs 


'^l 


A2 


^1 




n 
p 


P 




tllXLG 


1 


50 


40 


80 


1 


0.2 


2.0 


0.5 


2 


20 


2 


50 


40 


60 


1 


0.2 


2.0 


0.5 


2 


20 


3 


100 


80 


120 


1 


0.2 


2.0 


0.7 


2 


20 


4 


100 


90 


110 


1 


0.2 


2.0 


0.7 


2 


20 


5 


50 


40 


80 


1 


0.2 


1.5 


0.7 


2 


20 


6 


50 


40 


60 


1 


0.2 


1.5 


0.7 


2 


20 


7 


50 


45 


55 


1 


0.2 


2.0 


0.5 


2 


20 


8 


100 


95 


105 


1 


0.2 


2.0 


0.5 


2 


20 


9 


150 


140 


160 


1 


0.2 


2.0 


0.5 


2 


20 


10 


150 


100 


190 


1 


0.2 


2.0 


0.5 


2 
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of each experiment. In Table [T] "svrs" is the number of servers {nt), "alter" is the time length 
for which each arrival rate lasts, and "time" is the end time of our analysis. We already recog- 
nize that the method assuming measure zero works well when it does not linger too long near the 
non-differentiable points. For comparison, therefore, our experiments contain several cases where 
the system does linger relatively long around those points as well as the cases where it does not. 
Experiments 1-4 are intended to see the effects of "lingering" around non-differentiable points. We 
change f5t = 13 and pt = p as well as the arrival rates in experiments 5-8 to see the effects of other 
parameters. In fact, from the other experements not listed in Table [T] it turns out that changing 
other parameters does not affect estimation accuracy significantly. Experiments 9 and 10 are set 
to observe how larger arrival rates and number of servers affect estimation accuracy along with 
the lingering around the non-differentiable points by increasing both of them. Here we explain the 
overall results: for the details of numerical results, see Table |2]|6] in Appendix |B| Similar to the 
results in Figures [2] and [Sj we observe that lingering does debase the quality of approximations 
significantly when assuming measure zero. On the other hand, we see that our proposed method 
provides excellent accuracy in both mean and covariance values. Even if we increase both arrival 
rates and number of servers, we notice that lingering still affects estimation accuracy significantly 
when assuming measure zero but it does not in our proposed method. Figure |4] illustrates the 
average percentile difference of both methods against the simulation. The average difference are 
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Average difference (%) between each model and simulation 
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E[x1(t)] E[x2(t)] Var[x1(t)] Cov[x1 (t),x2(t)] Var[x2(t)] 

Performance measures 



Figure 4: Average difference against simulation 

obtained by averaging all differences in the tables, so it does not provide the absolute comparison 
between two methods. However, from Figure |4j we intuitively notice that our proposed method 
shows great effectiveness relative to the method assuming measure zero. 

6 Additional Applications 

We applied our proposed method to a wide variety of non-stationary state-dependent queueing 
systems. Since our proposed method is based on the adjusted fluid model, we observe that the 
mean queue lengths are accurate in all the systems. Also, when the rate functions are smooth or 
Gaussian density approximation is accurate or both, our adjusted diffusion model also provides 
accurate results. Due to space restrictions and our perception of how much value those cases 
would add, we have omitted presenting them here. Instead, we focus on scenarios with non- 
smooth rate functions where we conjecture Gaussian density would perhaps be inaccurate. We 
specifically consider two such applications not to showcase the effectiveness of our methodology, hut 
to illustrate that there is room for improvement for researchers in future to consider. We would 
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Figure 5: Multiclass priority queue with preemptive policy, (author?) [6] 

like to nonetheless point out that to the best of our knowledge our approximations are still more 
accurate than those in the present literature. In particular, we consider multiclass preemptive 



priority queues (Section |6.1[) and peer networks (Section 6.2). For these applications, we do not 



provide as much detail as the multiserver queues with abandonments and retrials in Section [5] and 
show results graphically. 

6.1 Multiclass preemptive priority queues 

In this section, we consider a multiclass priority queue (see Figure |5]) with preemptive policy which 
in fact is described in (author?) [6]. It is crucial to notice that (author?) [6^ does not numerically 
solve this example. However, we will not only use our methodology but also extend the method 
assuming measure zero in (author?) [7j for this case. 

Explaining the priority queue briefly, there are c number of classes of customers. The class i 
customers arrive to the system with rate XI and are served by available servers among rit number 
of servers with rate at time t. If a class i customer arrives and there is no available server, then 
the highest class customer (i.e. lowest priority) in service is pushed back to the queue and the class 
i customer is served. If there is no higher class customer, then the class i customer waits in the 
queue. 

In our numerical study, we use two classes of customers. Let X{t) = {xi{t),X2{t)^ be the state of 
the system at time t where xi{t) and X2{t) are the number of customers of class 1 and 2 respectively. 
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Then, X{t) is the solution to the following integral equations: 



We set the number of servers {nt) to be 200. The arrival rate of class 1 customers {\\) is alternating 
between 120 and 200 every two time units and the arrival rate of class 2 customers (A^) is 20. The 
service rates of both class customers are 1, i.e. = ^1 = 1. We conduct 5,000 simulation runs 
and obtain mean values by averaging them. 



n=200A,,=120,)i,2=200,A2=20,n=1,t=20, alter=2 



n=200,i„=120,i,2=200,)i2=20, n=l, t=20, alter=2 
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(a) Mean numbers by assuming measure zero (b) Mean numbers by our proposed method 

Figure 6: Comparison of mean values, E[X{t)\ 



Figures [6] and [7] show the comparison between the method assuming measure zero and our proposed 
method against simulation. We see that both methods work well for the mean value of xi{t). 
However, though not immediately obvious from Figure [6| there is 5-15% difference for the mean 
value of X2{t) when using the method assuming measure zero while our proposed method shows 
great accuracy. For the covariance matrix, our proposed method outperforms the method assuming 
measure zero as seen in Figure [Tj However, underestimation of variance of X2{t) is observed in our 
proposed method. Here, we explain our conjecture on the underestimation of variance. As described 
in Section|4| we utilize Gaussian density to obtain new rate functions, •)'s. In this example, we 
observe from our numerical experiments that empirical density is not close to Gaussian density when 
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(a) Covariance matrix by assuming measure zero (b) Covariance matrix by our proposed method 

Figure 7: Comparison of covariance matrix entries, , Cov^X(t), X{t)~^ 
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Figure 8: Peer networks 

the fluid limit stays near a non-difl'erentiable point. Our conjecture is that asymmetry of empirical 
density, unlike Gaussian density, causes larger values of covariance matrix entries. However, note 
that although it does affect the estimation of the covariance matrix (usually underestimation), it 
does not affect the estimation of mean values of the system state significantly, i.e. we still have the 
accurate estimation of mean values. 



6.2 Peer networks 

Figure [8] illustrates the queueing system we consider in this section. Motivated by peer networks 
with centralized controller (frequently encountered in multimedia content delivering industry), we 
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consider the following scenario of a queueing system. 

At time t, customers arrive to the system at rate At. Customers in the queue are sent to available 
active servers. The service rate of each server is fit- After a customer is fully served, s/he becomes 
a new active server. Each server serves customers for a random amount of time with mean l/9t, 
and then either becomes inactive with probability pt or leaves the system with probability 1 — pt. 
Each inactive server spends a random amount of time with mean and returns to be active. 
Note that only active servers can serve customers. 

Let X{t) = X2(i), X3(t)) be the state of the system at time t where xi{t), X2{t), and X3{t) are 

the number of customers, active servers, and inactive servers respectively. Then, X(t) is obtained 
by solving the following integral equation: 



where Yi{-), Y2{-), Y3(-), l4(-), and Y5{-) are independent Poisson processes with rate 1 correspond- 
ing to customer arrival, service completion, server's going up, server's leaving, and server's going 
down respectively. Similar settings are found in the previous research studies (? ], and ? ]). They, 
however, used Poisson processes (i.e. constant rate functions) to construct the system model and 
did not consider the time-varying properties. Furthermore, they mainly focused on the steady state 
analysis not providing in-depth analysis of the transient behavior of the system. Therefore, for this 
application, we can say that we provide more general settings than previous research in that we 
address time-varying rate functions and provide performance measures through the entire lifespan 
of the system. In this section, however, we provide the analysis of transient period when the system 
does not have enough servers to serve customers. This transient analysis is important since the 
issues on quality of service may arise during this period. We will provide details of the analysis of 
peer networks in a forthcoming paper. 
Figure lo] illustrates mean values of customers active servers (i? [2:2(0] )) inactive 
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(a) Covariance matrix by assuming measure zero (b) Covariance matrix by our proposed method 

Figure 10: Comparison of covariance matrix entries, Cov[^X(t), X(t)'j 
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servers (ii^[x3(t)]) over time. We apply both methods (i.e. the method assuming measure zero 
and our proposed method) and compare them with the simulation result. We use parameters with 
Xt = 400, = 2,9t = 0.3, 7t = 0.5, pt = 0.9 for t > as well as xi(0) = X3(0) = and X2(0) = 10. 
We conduct 5,000 simulation runs and obtain mean values by averaging them. In this numerical 
example, we want to see what happens when the fluid limit becomes close to the critically loaded 
time point. Figure |9] (a) shows comparison between simulation and the method assuming measure 
zero. We see that it works well when the fluid limits are far from the critically loaded time point. 
However, as it becomes closer to that point, the method assuming measure zero shows difference 
from the simulation result. On the other hand, when we apply our proposed method, as seen in 
Figure [9] (b), it provides the almost exact estimation even if the fluid limits are close to the critically 



loaded point. Figure 10 is the graph of the covariance matrix entries of the system over time. We 



observe sharp spikes in the method assuming measure zero (Figure 10 (a)). Note that the time 
when the extreme points of sharp spikes occur is exactly same as the time when the fluid limit 
hits the critically loaded time point in Figure [9] (a) which is the non-differentiable point. However, 
when we apply our proposed method, we have no sharp spikes at all and the covariance matrix 



entries are quite close to the simulation result as seen in Figure 10 (b). Thus we believe that our 
proposed method works well even under this complex scenario. 



7 Conclusion 

In this paper, we initially explain the fluid and diffusion models used in analysis of state-dependent 
queues and show potential problems that one faces in balancing accuracy and computational 
tractability. The first problem is from the fact that expectation of a function of a random vector X 
is not equal to the value of the function of the expectation of X. Therefore, unless they are equal or 
close, the fluid model may not provide an accurate estimation of mean values of the system state. 
The second problem is caused by non-differentiability of rate functions which prevents applying 
the diffusion model in (author?) [5j. Therefore, addressing these problems is quite important in 
order to develop accurate approximations as well as to achieve computational feasibility. For that, 
we proposed a methodology to obtain the exact estimation of mean values of system states and an 
algorithm to achieve computational tractability. 



28 



The basic idea of our approach is to construct a new stochastic process which has the fluid hmit 
exactly same as the mean value of the system state. We proved that if rate functions in the original 
model satisfy the conditions to apply the fluid model, rate functions in the constructed model also 
satisfy those conditions. Therefore, we can apply the adjusted fluid model if we can apply the 
existing fluid model. It turns out that there is, in general, no computational method to obtain the 
adjusted fluid model exactly and hence we utilize Gaussian density to approximate it. By using 
Gaussian density, we see that rate functions in the constructed model are smooth and we are able 
to apply the diffusion model in (author?) [5j even if we could not apply it to the original process. 
To validate our proposed method, we provide several numerical examples of non-stationary state- 
dependent queueing systems. In the examples, we observe that our proposed method shows great 
accuracy compared with the method assuming measure zero (which is the only other way in the 
literature, to the best of our knowledge, that provides computational tractability). Due to space 
restriction, we have not shown all examples where our method works well. We, however, ob- 
serve that when the Gaussian density assumption is inaccurate, especially near non-smooth points, 
our methodology needs further investigation for the covariance matrix. We conjecture that this 
phenomenon is from the gap between empirical and Gaussian density. To address this, one can 
investigate the properties of speciflc rate functions that affect the shape of empirical density and 
can devise a new algorithm flnding (jrj(-, •)'s from other density functions in the future. 
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Appendices 



A Derivation of gi{t,x)''s 

For fixed to > 0, let n = ntg, m = fij^, (3 = (3to, P = Pto, xi = xi(to) ~ N{zi,al), and X2 = X2(to) 
N{z2,a2)- For z = (21,^2)', we have 



g?.{toiz) = A n)] = /ii|£;[xila;i<n] + nPr[xi > n]| 

(x - Zif^ 



exp 



27r(Ti 

0"! X — Zl 
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+ nPr[xi > n] 



'27r J -00 

o 1 



'1 



— exp 
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(x - ziY' 



2af 



dx + ziPr[xi < n] + nPr[xi > n] 



2-irai 



{n - zif 
2al 



+ (zi — n)Pr{xi <n)+n 



Therefore, by making to > arbitrary, we have g^{t,x). 

Note gi{-, •) and g^{-, •) are same except a constant part with respect to x. Therefore, it is enough 
to derive 55(-, •)• We can show that 



55(^0, ^) = E[l3p[xi-n)+]= l5p{E[xi'^n\-n] 
= /3p|£'[xil^j>„] + nPr[xi < n] - n| 

(x - zi)2 



f3p 
I5p 
(3p 



n \2lTai 

00 



exp 



2af 



dx + nPr[xi < n] — n 



-cri 
/2^ 



n (^1 



exp 



X — Zi 
2 ^^P 



2 — ) '^^ + ziPr[xi > n] + nPr[xi < n] — n 



2cr 



27r(Ti 



(ra - ziY' 

2g\ 



+ (zi — n)Pr{x\ > n) 



Therefore, by making to > arbitrary, we have gf5(t, x). 



B Numerical results for Section [5] 
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Table 2: Estimation of over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


() 




8 


9 


10 


il 


12 


i:i 


11 


15 


1 


proposed 


6.52 


0.98 


-3.39 


-1.07 


-3.05 


-0.40 


0.91 


0.25 


-0.69 


-0.01 


rneas. 


4.42 


0.82 


-3.63 


-1.94 


-3.60 


-0.23 


0.75 


-0.15 


-2.59 


0.11 


2 


proposed 


2.69 


0.44 


-3.13 


-0.82 


-1.08 


-0.32 


0.48 


0.15 


-0.46 


-0.05 


meas. 


3.35 


-0.42 


-2.92 


-1.64 


-1.18 


-1.01 


0.85 


-0.44 


-0.36 


-0.60 


3 


proposed 


2.33 


0.28 


-3.11 


-1.01 


-1.36 


-0.39 


0.10 


-0.02 


-1.68 


-0.15 


rneas. 


2.34 


-0.42 


-2.67 


-1.55 


-1.49 


-1.00 


0.52 


-0.49 


-0.54 


-0.53 


4 


proposed 


1.18 


0.14 


-1.54 


-0.30 


-0.01 


0.12 


0.22 


0.22 


-0.10 


-0.02 


rneas. 


0.65 


-0.96 


-1.98 


-1.32 


-0.94 


-0.95 


0.04 


-0.64 


-0.61 


-0.94 


5 


proposed 


7.04 


1.36 


-3.67 


-0.69 


-1.38 


-0.57 


0.80 


0.23 


-2.82 


-0.63 


meas. 


5.55 


1.04 


-3.20 


-0.93 


-1.31 


-0.53 


0.46 


0.06 


-1.22 


-0.18 


6 


proposed 


3.61 


0.76 


-3.05 


-1.13 


-0.67 


0.18 


1.12 


0.20 


-0.95 


-0.25 


meas. 


2.53 


-0.07 


-3.01 


-1.72 


-1.46 


-0.43 


0.60 


-0.47 


-1.57 


-0.80 


7 


proposed 


1.93 


0.65 


-1.06 


-0.25 


-0.63 


0.17 


0.12 


-0.21 


-0.65 


-0.20 


meas. 


0.50 


-0.86 


-2.07 


-1.51 


-1.04 


-0.73 


-0.47 


-1.07 


-0.63 


-0.76 


8 


proposed 


0.72 


0.07 


-0.46 


0.04 


-0.04 


-0.14 


0.42 


-0.07 


-0.48 


-0.01 


meas. 


0.04 


-0.98 


-1.40 


-0.91 


-0.57 


-0.85 


-0.13 


-0.69 


-0.73 


-0.46 


9 


proposed 


0.81 


0.25 


-0.96 


-0.25 


-0.11 


-0.09 


0.38 


-0.06 


-0.24 


-0.02 


meas. 


0.53 


-0.50 


-1.31 


-0.88 


-0.34 


-0.61 


0.17 


-0.51 


-0.06 


-0.32 


10 


proposed 


6.44 


1.18 


-4.73 


-1.73 


-2.21 


-0.45 


0.30 


-0.01 


-1.10 


-0.11 


meas. 


6.46 


0.77 


-3.83 


-1.62 


-2.84 


-0.83 


0.84 


0.00 


-2.77 


-0.60 



Table 3: Estimation of i?[x2(t)] over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


-2.00 


3.50 


2.36 


-0.53 


0.57 


-1.00 


-0.99 


-0.30 


-0.44 


-0.76 


meas. 


11.68 


12.60 


7.38 


5.88 


11.64 


8.18 


5.24 


6.29 


10.47 


7.82 


2 


proposed 


-2.22 


2.71 


1.90 


-2.44 


-0.94 


-1.82 


-0.91 


-0.10 


-0.38 


-0.76 


meas. 


45.00 


53.07 


33.49 


37.12 


41.73 


44.51 


31.55 


37.48 


40.57 


43.21 


3 


proposed 


-2.49 


1.88 


1.00 


-3.58 


-2.09 


-3.32 


-3.08 


-3.01 


-3.15 


-4.02 


meas. 


28.64 


37.65 


19.44 


21.88 


24.73 


28.38 


16.37 


21.45 


22.77 


26.96 


4 


proposed 


0.24 


2.66 


1.35 


-1.68 


-0.91 


-0.19 


0.25 


0.45 


0.02 


-0.53 


meas. 


67.95 


69.81 


47.03 


51.81 


56.69 


59.66 


45.16 


50.75 


54.48 


57.27 


5 


proposed 


-1.01 


4.41 


3.16 


-0.05 


1.55 


0.57 


-0.58 


-0.07 


-0.09 


-1.63 


meas. 


9.61 


12.28 


7.50 


5.73 


11.28 


9.42 


5.36 


6.00 


9.51 


8.02 


6 


proposed 


-2.63 


2.48 


2.23 


-2.39 


-1.32 


-0.84 


-0.12 


1.04 


0.89 


-0.04 


meas. 


44.23 


51.84 


32.45 


35.11 


39.27 


43.72 


31.00 


35.94 


38.83 


41.83 


7 


proposed 


0.33 


3.42 


3.00 


1.01 


0.70 


0.25 


0.41 


0.71 


0.27 


-0.17 


meas. 


78.08 


78.96 


60.84 


64.86 


69.59 


71.19 


59.15 


63.20 


67.42 


68.95 


8 


proposed 


2.81 


3.03 


2.40 


1.45 


1.29 


0.58 


-0.08 


0.41 


0.12 


-1.11 


meas. 


92.68 


90.60 


73.97 


77.06 


80.98 


81.48 


70.82 


74.24 


77.96 


78.55 


9 


proposed 


-0.86 


1.25 


1.44 


-0.77 


-0.19 


-0.18 


0.08 


0.84 


0.42 


0.41 


meas. 


80.15 


79.90 


57.59 


62.03 


67.09 


68.91 


55.09 


59.98 


64.19 


66.14 


10 


proposed 


-2.07 


().()2 


:i.79 


-2.r)0 


-0. 19 


-2.18 


-1.99 


-l.:!5 


-l.:.i8 


-1.21 


meas. 


8.53 


23.91 


10.73 


8.77 


10.78 


13.05 


5.67 


8.96 


9.13 


11.69 
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Table 4: Estimation of Far [a;i(t)] over time; difference from simulation 



ExporiiiK^nts 


Time (t) 


# 


mot hod 


(i 




8 


9 


iO 


il 


12 


i;i 


11 


15 


1 


proposed 


6.94 


0.94 


-1.92 


-2.02 


-3.66 


-0.20 


1.70 


-1.02 


0.49 


2.89 


mcas. 


-11.03 


2.93 


-1.93 


17.31 


-24.44 


1.42 


1.66 


14.89 


-19.73 


4.16 


2 


proposed 


2.84 


3.83 


-6.05 


-0.10 


-0.50 


4.24 


1.62 


2.67 


-1.02 


1.62 


meas. 


-6.28 


16.69 


6.76 


-14.45 


-12.97 


17.62 


12.90 


-11.61 


-14.51 


15.29 


3 


proposed 


4.15 


2.09 


-0.60 


2.15 


-6.57 


0.50 


-3.10 


1.15 


2.76 


3.74 


meas. 


-0.56 


13.38 


7.30 


-8.74 


-12.97 


11.92 


4.53 


-10.16 


-2.13 


14.22 


4 


proposed 


-0.52 


-4.36 


-2.81 


3.07 


-0.03 


2.96 


0.79 


1.27 


3.30 


0.35 


meas. 


-16.38 


11.18 


14.13 


-12.86 


-17.81 


17.93 


16.94 


-15.33 


-14.05 


15.32 


5 


proposed 


6.83 


-0.22 


-2.49 


0.09 


-1.67 


-3.27 


1.71 


-4.14 


-0.55 


1.98 


meas. 


-2.30 


1.03 


-1.69 


10.07 


-10.59 


-1.97 


1.43 


5.09 


-7.69 


3.42 


6 


proposed 


5.22 


0.62 


-6.25 


-0.81 


-4.32 


-1.95 


4.41 


1.97 


-0.61 


4.93 


meas. 


-1.19 


7.72 


1.39 


-7.42 


-11.70 


5.61 


10.28 


-4.33 


-7.73 


12.15 


7 


proposed 


2.91 


-2.29 


-1.04 


0.92 


0.21 


0.18 


3.14 


-1.10 


4.36 


2.28 


meas. 


-17.83 


14.52 


18.27 


-16.55 


-22.37 


17.07 


20.88 


-18.77 


-18.14 


19.01 


8 


proposed 


-1.79 


0.65 


-0.43 


0.83 


3.35 


-0.71 


3.63 


2.10 


1.85 


0.72 


meas. 


-26.38 


16.44 


21.26 


-18.37 


-22.66 


16.73 


23.72 


-17.42 


-25.80 


18.25 


9 


proposed 


0.62 


-0.86 


-0.83 


3.53 


3.36 


5.09 


1.52 


1.71 


2.73 


-1.37 


meas. 


-17.84 


13.72 


17.09 


-14.12 


-17.40 


19.78 


18.07 


-16.57 


-19.03 


14.36 


10 


proposed 


4.48 


-0.32 


-9.84 


1.26 


-4.24 


3.37 


1.32 


1.00 


0.22 


1.15 


meas. 


4.12 


7.27 


-6.22 


-2.69 


-5.55 


10.87 


3.68 


-3.27 


-2.12 


8.98 



Table 5: Estimation of Cov[xi{t),X2{t)~\ over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


-3.03 


-3.27 


4.75 


3.10 


-1.72 


-3.63 


0.39 


-4.00 


-4.15 


0.91 


meas. 


25.05 


-3.88 


4.43 


-6.75 


15.78 


-4.69 


-0.30 


-11.23 


4.26 


-2.03 


2 


proposed 


-6.76 


7.23 


2.87 


-4.73 


11.50 


-0.47 


6.09 


5.48 


5.86 


4.82 


meas. 


29.60 


-6.12 


-6.56 


3.81 


36.90 


-13.05 


-1.41 


12.63 


31.69 


-5.53 


3 


proposed 


-6.74 


-2.24 


4.53 


-9.51 


-28.97 


-3.44 


-2.57 


-6.43 


5.52 


-3.76 


meas. 


25.67 


-15.55 


0.47 


6.26 


6.46 


-15.26 


-6.35 


8.68 


31.03 


-13.44 


4 


proposed 


-0.01 


-13.57 


-8.53 


-12.42 


-9.73 


-0.00 


-10.28 


-16.74 


-1.43 


-11.64 


meas. 


58.61 


-29.39 


-21.29 


10.14 


44.87 


-14.34 


-22.28 


5.51 


46.98 


-26.05 


5 


proposed 


-7.19 


-0.18 


0.13 


2.88 


7.13 


-8.20 


-0.42 


1.07 


4.17 


-1.88 


meas. 


19.73 


-4.03 


-1.03 


-12.97 


26.24 


-11.11 


-1.32 


-12.31 


20.40 


-4.25 


6 


proposed 


2.91 


4.90 


2.63 


-7.64 


-7.99 


1.17 


2.14 


3.80 


6.94 


7.92 


meas. 


37.93 


-15.54 


-15.96 


-3.32 


25.88 


-18.86 


-15.00 


6.43 


34.79 


-10.85 


7 


proposed 


-4.38 


-1.88 


0.97 


-14.36 


3.15 


-1.95 


-1.19 


-1.08 


-4.77 


-0.46 


meas. 


52.91 


-16.88 


-16.53 


-3.19 


43.26 


-15.54 


-16.02 


6.73 


34.99 


-12.29 


8 


proposed 


-20.99 


-4.21 


-6.33 


-4.51 


3.12 


-3.43 


-1.85 


-6.78 


-1.81 


-0.79 


meas. 


64.94 


-8.85 


-22.74 


13.30 


51.79 


-11.89 


-15.66 


7.80 


44.16 


-8.72 


9 


proposed 


-15.01 


-6.15 


-6.27 


3.33 


-0.25 


2.45 


-6.00 


-7.57 


-6.93 


-6.74 


meas. 


55.84 


-12.34 


-17.97 


19.55 


45.76 


-4.17 


-15.42 


7.67 


37.97 


-12.70 


10 


proposed 


-18.70 


7.57 


-:i.70 


-1.76 


8.01) 


-(). 13 


-2.8(i 


-0.03 


2.9:) 


-1.11 


meas. 


-21.43 


-2.63 


-5.67 


-0.67 


8.99 


-15.71 


-4.40 


3.66 


4.18 


-9.87 
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Table 6: Estimation of Far [a;2(i)] over time; difference from simulation 



Experiments 


Time (t) 


# 


type 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


1 


proposed 


-2.15 


3.52 


1.48 


-0.72 


-0.34 


-0.45 


0.78 


1.59 


1.31 


0.83 


meas. 


6.74 


5.31 


2.34 


-8.01 


5.46 


1.00 


1.84 


-2.78 


2.81 


-1.20 


2 


proposed 


1.29 


9.81 


8.50 


3.31 


7.72 


6.70 


6.05 


5.84 


5.42 


5.91 


rncas. 


7.06 


14.88 


-6.56 


-0.09 


17.07 


12.12 


-3.90 


5.44 


16.94 


13.19 


3 


proposed 


-5.60 


2.13 


-0.71 


-5.14 


-1.93 


-3.21 


-2.71 


-1.18 


-0.61 


-0.79 


meas. 


-2.22 


1.5.'i 


-10. ■'57 


-5.1 1 


i.Ol 


1.00 


-8.89 


().()() 


(). 15 


5.96 


4 


proposed 


5.71 


8.34 


2.63 


-2.01 


-0.83 


2.01 


-0.03 


-0.27 


0.51 


2.66 


meas. 


28.49 


26.24 


-13.87 


-3.06 


15.93 


15.77 


-12.78 


0.50 


17.18 


17.62 


5 


proposed 


-0.97 


4.30 


1.25 


-0.07 


3.22 


3.50 


-0.22 


-0.16 


1.50 


0.76 


meas. 


3.67 


5.28 


1.60 


-4.35 


7.92 


6.18 


1.65 


-2.36 


5.70 


4.03 


6 


proposed 


2.23 


11.10 


8.33 


2.94 


4.21 


4.30 


1.38 


2.79 


3.63 


5.03 


meas. 


11.13 


21.61 


-3.19 


-0.25 


12.83 


14.35 


-6.15 


1.31 


12.95 


14.75 


7 


proposed 


5.23 


7.48 


5.57 


1.60 


2.20 


4.24 


3.45 


4.88 


5.03 


4.33 


meas. 


33.03 


27.56 


-16.08 


-4.39 


21.67 


19.11 


-11.03 


2.36 


24.85 


19.64 


8 


proposed 


10.11 


7.03 


3.99 


2.44 


3.47 


2.45 


2.53 


2.25 


1.73 


2.30 


meas. 


62.03 


46.52 


-13.63 


0.58 


30.10 


23.25 


-13.94 


-0.09 


26.38 


20.60 


9 


proposed 


8.18 


7.49 


3.22 


0.53 


3.83 


4.55 


3.14 


4.24 


4.90 


5.28 


meas. 


39.93 


31.88 


-18.20 


-4.36 


22.39 


18.66 


-12.73 


1.21 


23.00 


18.98 


10 


proposed 


-0.34 


12.31 


5.05 


-2.01 


1.38 


1.13 


-3.01 


-3.64 


-4.35 


-1.66 


meas. 


-5.73 


7.15 


-1.91 


-3.68 


0.93 


-2.52 


-8.00 


-4.34 


-3.94 


-5.72 
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